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An exact expression for the leading-order (LO) gluon distribution function G(x,Q 2 ) = xg(x,Q 2 ) 
from the DGLAP evolution equation for the proton structure function F^ v (x,Q 2 ) for deep inelastic 
7*p scattering has recently been obtained [M. M. Block, L. Durand and D. W. McKay, Phys. 
Rev. D79, 014031, (2009)] for massless quarks, using Laplace transformation techniques. Here, we 
develop a fast and accurate numerical inverse Laplace transformation algorithm, required to invert 
the Laplace transforms needed to evaluate G(x, Q 2 ), and compare it to the exact solution. We obtain 
accuracies of less than 1 part in 1000 over the entire x and Q 2 spectrum. Since no analytic Laplace 
inversion is possible for next-to-leading order (NLO) and higher orders, this numerical algorithm 
will enable one to obtain accurate NLO (and NNLO) gluon distributions, using only experimental 
measurements of FJ P (x,Q 2 ). 

(N ■ 

I. INTRODUCTION 

| The quark and gluon distributions in hadrons play a key role in our understanding of Standard Model processes, 
in our predictions for such processes at accelerators, and in our searches for new physics. In particular, accurate 
knowledge of gluon distribution functions at small Bjorkcn x will play a vital role in estimating backgrounds, and 
hence, our ability to search for new physics at the Large Hadron Collider. 

The gluon and quark distribution functions have traditionally been determined simultaneously by fitting experi- 
Oh| mental data (mainly at small x) on the proton structure function F^ix^Q 2 ) measured in deep inelastic ep (or 7*p) 
scattering, over a large domain of values of x and Q 2 . The process starts with an initial Qq, typically in the 1 to 2 
GeV 2 range, and individual quark and gluon trial distributions parameterized as functions of x. The distributions are 
evolved to larger Q 2 using the coupled integral-differential DGLAP equations P, d, UJ], and the results used to predict 
the measured quantities. The final distributions are then determined by adjusting the input parameters to obtain a 
best fit to the data. For recent determinations of the gluon and quark distributions, see 0, [fill, Q- 

This procedure is rather indirect, especially so in the case of the gluon: the gluon distribution G(x,Q 2 ) does not 
appear in the experimentally accessible quantity F^ p (x, Q 2 ), and is determined only through the quark distributions 

■ in conjunction with the evolution equations. It is further not clear without detailed analysis [H, HI, 0, [13] how sensitive 
. ' | the results are to the parameterizations of the initial parton distributions, or how well the gluon distribution is actually 

• determined. 

' In a recent paper, Block, Durand and McKay (BDM) [ll[ used a Laplace transformation technique to obtain 
C^) . a leading order (LO) analytic gluon distribution function G(x,Q 2 ) = xg(x,Q 2 ) for massless quarks, directly from 
a global parameterization of the data on F2 P (x,Q 2 ). The method uses only the LO DGLAP evolution equation 
El for F2 P {x, Q 2 ) in the usual approximation in which the active quarks are treated as massless. In contrast 

■ to previous methods for determining G(x 7 Q 2 ), it does not require knowledge of the separate quark distributions in 
the region in which structure function data exist, nor does it require the use of the evolution equation for G(x,Q 2 ), 
both considerable simplifications. In essence, the authors transform the LO DGLAP equation from Bjorken x-space 
to w-space, where v = ln(l/x), and then Laplace transform the resulting equation. After solving for the Laplace 
transform g(s,Q 2 ), they are able to analytically invert the Laplace transform back into v-space, and eventually, to 
x-space. Unfortunately, because of the considerable complexities of next-to-leading order (NLO) splitting functions 
needed in the NLO DGLAP evolution equation [l], 0, 9 for F^ p (a;, Q 2 ) , the method can not be extended to NLO 
gluon distributions, because of the impossibility of analytically inverting the required Laplace transform. 

The purpose of this note is to derive a new and fast algorithm for accurate numerical inversion of Laplace transforms, 
so that the BDM method [ll[ for direct evaluation of gluon distributions from knowledge of F^ix, Q 2 ) can be extended 
to NLO (and NNLO) easily and accurately. Mathematically, Laplace inversion is an "ill-posed" problem and great 
care must be taken to insure its reliability for arbitrary Laplace transforms 0. In this case, we can check our 
numerical Laplace inversion routine directly by comparing its results to the exact LO gluon solution of Ref . [Til ] . 
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II. THE EXACT LO SOLUTION 



The LO DGLAP equation for the evolution of the proton structure function K^ p (:r,<3 2 ) for 4 massless quar 
(u, d, s, c) can be written as 



dF^{x,Q 2 
d\nQ 2 



47T 



dz 



(z,Q 2 )K qq (^) 
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dz 



G{z,Q 2 )K q 



where K qq (x) and K„„ (x) are the LO splitting functions and a s is the strong running coupling constant. 
Following BDM [nf, we introduce F% p (x, Q 2 ) 



dF^jx, Q 2 
d\nQ 2 



47T 



1 dz 



71 r2 



\z,Q 2 )K qq {^). 



We finally write the DGLAP equation for the evolution of F*[ (x,Q ) as 



where 



x / G{z,Q 2 )K qg (-\ — = F{x,Q% 



and the LO g — > q splitting function is given by 

K qg (x) = l-2x + 2x 2 . 
At this point, BDM introduce the coordinate transformation 

v = ln(l/x), 

and define functions G, K qg , and T in u-space by 

G(v,Q 2 ) = G{e-\Q 2 ) 

K qg {v) ee K qg (e- V ) 
Hv.Q 2 ) = F(e~ v ,Q 2 ). 



Explicitly, from Eq. ([5]), we see that 



Since 



K qg (v) = 1 - 2e~ v + 2e 



T{v,Q 2 ) = / G(w ) Q 2 )e- ( - v - w ^K qg (v-w)dw, 
Jo 

BDM note that in i>-space, the DGLAP equation of Eq. ([3]) can be written as the convolution integral 



T{v,Q 2 



G(w,Q 2 )H(v-w) dw, 



where 



H(v) ee e~ v K qg {v) 

= e~ v - 2e- 2v + 2e~ 3v . 

The Laplace transform of H(v) is given by h(s), where 

/>oo 

h(s) ee C[H{v); s}= H(v)e~ sv dv, with H{v) = 0, v < 0. 
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The convolution theorem for Laplace transforms allows us to rewrite the right hand side of Eq. (|10[) as a product of 
their Laplace transforms g(s) and h(s), so that the Laplace transform of Eq. (|10[) is given by the algebraic equation 



f(s,Q 2 ) = g(s,Q 2 )xh(s). 
Solving Eq. (|13p for g, the Laplace transform of the gluon distribution function in s-space is given by 

g(s,Q 2 ) = f(s,Q 2 )/h(s). 



(13) 



(14) 



In general, one is not able to calculate the inverse transform of g(s, Q 2 ) explicitly, if only because f(s, Q 2 ) is determined 
by a numerical integral of the experimentally-determined function T(v,Q 2 ). However, regarding g(s,Q 2 ) as the 
product of the two functions f(s,Q 2 )a,nd h^ 1 (s) and taking the inverse Laplace transform using the inverse of the 
convolution theorem and the known inverse £~ 1 [f(s, Q 2 ); v] = P(v, Q 2 ), BDM find the analytic solution 



G(v,Q 2 ) = C-'ifis^^xh-^s^v] 



(15) 



The calculation of h(s) and the inverse Laplace transform of h 1 (s) are straightforward, and the answer for LO, 
albeit singular, is given by BDM in terms of the Dirac delta function as 



£-* [ft" 1 («);«] = 36(v)+5'(v)-er 3v / 2 
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Thus, again using the convolution theorem, BDM find that the analytic solution for the LO gluon distribution G(v, Q 2 ) 
for massless quarks is 



G{v,Q 2 ) = 3F(v,Q 



2 | df(v,Q 2 ) 
dv 

T( Wl Q 2 )e-^ v - 



sm 



V7, 

— {v-w) 



2 cos 



V7, 

— {v-w) 



dw. 



(16) 



The BDM solution depended critically upon the ability to find the analytic inverse Laplace transform of ft,~ 1 (s), which 
was possible for the LO case of the splitting function K qg (x). 

Unfortunately, for NLO or higher, the splitting function is so complicated that an analytic inversion for the NLO 
h^ 1 is impossible — see Floratos et al. [13| for the NLO MS splitting function that is required. For this case, we 
must be able to find a numerical inversion of the Laplace transform g(s) of the equivalent of Eq. Q14p in order to 
obtain G(v, Q 2 ) and, ultimately, G(x, Q 2 ). The goal of this communication is to develop a suitable numerical Laplace 
inversion algorithm. 



III. NUMERICAL INVERSION OF LAPLACE TRANSFORMS 

In order to simplify our notation, we will now suppress the explicit dependence of G{v,Q 2 ) on Q 2 , writing it as 
G{v). If the Laplace transform of g(s) = C[G(v); s] = J °° G{v)e~ vs dv, where G{v) = for v < 0, then the inverse 
Laplace transform is given by the complex Bromwich integral 

i p c-\-i oo 

G{v) = C- 1 [g{s)-M = i r -. / 9(sV s ds, (17) 

where the real constant c is to the right of all singularities of <?(s). We will further assume that we have made an 
appropriate coordinate translation in s so that c = 0, so that the equation is written as 

i r +ioo 

G(v) = C-'igis); v] = — / g(s)e v ' ds. (18) 

Our goal is to numerically solve Eq. (|L8|) . The inverse Laplace transform is essentially determined by the behavior of 
g(s) near its singularities, and thus is an ill-conditioned or ill-posed numerical problem. We suggest in this note a 
new algorithm that takes advantage of very fast, arbitrarily high precision complex number arithmetic that is possible 
today in programs like Mathematica [l4| . making the inversion problem numerically tractable. 
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First, we introduce a new complex variable z = vs and rewrite Eq. (|18p as 



1 r +i °° rz\ 

G ( v ) = 7^— / 9(-)e z dz. (19) 



2niv _ 

We next make a rational approximation to e z , using the partial fraction expansion 



e-«£— ' (20) 

^ Z-Oii 

which can be shown to have the following properties: 

1. Rc Qi > 0, so that its poles are all in the right-hand half of the complex plane. 

2. iV distinct complex conjugate pairs of the complex numbers (u>i, ai), such that the sum of the fc th pair. 



(21) 



is real for all real 



3. The expansion is identical to the Pade approximant with numerator equal to 27V — 1 and denominator equal to 
27V. 

4. The integrand vanishes faster than l/R as R — ► oo on the semi-circle of radius R that encloses the right hand 
half of the complex plane, since the approximation vanishes as l/R and g(s) that corresponds to a non-singular 
G(v) also vanishes for R — > oo. 

Since g(z/v) must vanish for z — > oo and our approximation for e z in Eq. (|20[) vanishes for z — > oo, we can form a 
closed contour C by completing our integration path of the modified Bromwich integral in Eq. (| 1Q|) with an infinite 
half circle to the right half of the complex plane. As mentioned earlier, g(z/v) has no singularities in this half of the 
complex plane. It is important to note that this contour is a clockwise path around the poles of Eq. (|20p . which come 
from our approximation to e z . What we need is the negative of it, i.e., the contour — C which is counterclockwise, so 
that the poles are to our left as we traverse the contour — C. Therefore, we rewrite Eq. (fl9|) as 



1 f (z 



IN 



i=l 

2N 



Z ■ 



9 (-) — dz 



2ttw ' T_n 

2 — 1 ^ 

2 N 

= — VRe[ WiJ (ai/«)]. (22) 
u * — ' 

i— 1 

To obtain Eq. (|22p. the final approximation to G(v), we used Cauchy's theorem to equate the closed contour integral 
around the path — C to 2iri times the sum of the (complex) residues of the poles. Since the contour —C restricts 
us to the right-hand half of the complex plane, no poles of g(z/v) were enclosed, but only the 2iV poles a% of the 
approximation of e z . Using the properties cited above of the complex conjugate pairs — {u>i,cti) and (a>j,<5i) — after 
taking only their real part and multiplying by 2, we have simultaneously insured that G(v) is real , yet only have had 
to sum over half of the residues. 

Equation (|22[) has some very interesting properties: 

1. The 47V coefficients (a^Wj) arc complex constants that arc independent of v, only depending on the value of 
2N used for the approximation, so that for a given 27V, they only have to be evaluated once — in essence, they 
can be tabulated and stored for later use. 

2. The Laplace transform of v n is given by nl/s n+1 . where n is integer. Inserting G(v) = v n and g(s) = n\j ' s n+1 ^ 
into Eq. (|2"2"|) , we see that we have a set of 47V equations, 

N { \ 

-n!^2Ref-^ T J =1, n = 0, 1, . . . , 47V - 1, (23) 
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which also uniquely determine the 47V complex constants (a,, uii), i = 1,2,..., 27V, although evaluating them 
directly from Eq. (|23[) is virtually impossible numerically, considering the ill-posed nature of these equations. 

The true power of Eq. (|23|) that it shows that the inversion algorithm of Eq. ([22)) is exact when G(v) is a 
polynomial of order < 47V — 1, even though only N complex terms have to be evaluated in Eq. (|22ll . This is 
reminiscent of the situation using Gauss-Legendre integration of order N, where there are 27V constants, N 
Legendre zeroes and N weights, and the integration approximation is exact if the integrand is a polynomial of 
order < 27V- 1. 

3. The real parts of the residues u>i alternate in sign and are exceedingly large — even for relatively modest 27V, 
making round-off a potentially serious problem. Thus, exceedingly high precision complex arithmetic is called 
for, often requiring 60 or more digits. However, this is not a serious problem — either in speed or complexity of 
execution — for an algorithm written in a program such as Mathematica fl4| . 

A concise inversion algorithm in Mathematica that implements Eq. (|22[) is given in Appendix [X] and a one line 
Mathematica algorithm for finding a numerical Laplace transform is given in Appendix [B"l 

We will now test the accuracy of our numerical Laplace inversion algorithm by comparing its results with Eq. . 
the exact LO G(v, Q 2 ) of BDM. 

IV. COMPARISON OF EXACT SOLUTION AND NUMERICAL LAPLACE INVERSION RESULTS 

Using the Berger, Block and Tan [l5| fit to the experimental ZEUS [l6| data shown in Appendix [Cl we have 
calculated both the exact solution for G(v,Q 2 ) from Eq. (fT6|) and the completely numerical approximation for G(v) 
given in Eq. (|22[) . using our inverse Laplace transformation algorithm given in Appendix[A] For this purpose, we used 
2N = 8 and prec = 80 in the algorithm. The results for Q 2 = 5 GeV 2 and 100 GeV 2 are shown in Fig. Q] and Fig. 
[3J respectively. The blue points are the exact solution and the red curves are the numerical solution that uses our 
algorithm for the numerical inversion of Laplace transforms. As seen in both Figures, the agreement is striking over 
the entire v range, which corresponds to the x interval 5 x 10~ 7 < x < 1. At the highest v values (where numerical 
Laplace inversion approximations generally have the greatest problem), we find an accuracy of ~ 1 part in 5000, for 
Q 2 = 100 GeV 2 . This error reflects both the numerical inaccuracy associated with the numerical integration term in 
the exact solution G(v, Q 2 ) of Eq. (fTT)]) . as well as the inaccuracies associated with the numerical Laplace inversion 
routine of Appendix [X] and the numerical Laplace transformation routine of Appendix [Bj 

Another independent method of checking the numerical accuracv of the entire procedure is to go back to the original 
DGLAP equation from which we started, Eq. ([3]), i.e., 

x f G(z, Q 2 )K qg (-) ^ = F(x, Q 2 ), (24) 

J x \ Z / Z 

and numerically integrate its l.h.s., which depends on our numerical solution for G(v,Q 2 ) through G(x,Q 2 ) = 
G(\n(l/x),Q 2 ). We then compare it with the r.h.s., !F(x,Q 2 ), which is independently known for all x and Q 2 . 
This check becomes of primary importance when one does not have an analytical solution — the typical situation. In 
this case, it validated our conclusions about the numerical accuracy of our inversion procedure over the entire x and 
Q 2 domain. In particular, for Q 2 = 100 GeV 2 , the ratio of the l.h.s. to the r.h.s. of Eq. (|24[) is unity to 1 part in 
5000 in the x range 3 x 10~ 4 <i<3x 10 -2 , and never rises to more than ~ 1 part in 200 outside this range; these 
excursions from unity include the errors generated in the numerical integration of the l.h.s of Eq. (|24p . 

In contrast to the highly singular behavior of the analytic Laplace inverse of inversion of /i -1 (s) in Eq. ([To]) , the 
overall behavior of G(v,Q 2 ), the Laplace inversion of f(s,Q 2 )/h(s) in Eq. (|16p . is very smooth and therefore lends 
itself readily to numerical inversion, as seen by the excellent agreement shown in Fig. [T]and Fig. [5] between the exact 
solutions and the numerical inverse Laplace transforms. 

Although our inversion routine was specifically developed to invert Laplace transforms of gluon distributions, 
it clearly has a wide variety of applications in solving both integral and differential equations, which will not be 
commented on further in this communication. 



V. CONCLUSIONS 

We have achieved high numerical accuracy in obtaining LO gluon distributions from fits to experimental data, 
using a numerical Laplace inversion routine developed for this purpose, and have compared it successfully to the exact 
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FIG. 1: LO gluon distribution functions G(v,Q 2 ) vs. v = ln(l/x), for virtuality Q 2 = 5 GeV 2 . The blue points are the exact 
LO solution of Eq. (|16[) . The red curve is our numerical Laplace inversion solution. The agreement is excellent over the entire 
v range, which corresponds to the a>range, 5 x 10~ 7 < x < 1. 
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FIG. 2: LO gluon distribution functions G(v,Q 2 ) vs. v = ln(l/x), for virtuality Q 2 = 100 GeV 2 . The blue points are the exact 
LO solution of Eq. (|16[) . The red curve is our numerical Laplace inversion solution. The agreement is excellent over the entire 
v range, which corresponds to the rc-range, 5 x 10~ 7 < x < 1. 

analytic solution using the same fit. Since NLO and NNLO DGLAP solutions are not amenable to having analytic 
solutions — one can not analytically invert their h~ 1 (s) — this technique will allow us to find very accurate numerical 
gluon solutions directly from experimental data, without having first to find quark distributions by solving coupled 
DGLAP equations. Further, we will be able to verify the numerical accuracy of our solutions by direct substitution 
into their generating equations. 
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APPENDIX A: MATHEMATICA LAPLACE INVERSION ALGORITHM 



NInverseLaplaceTransf ormBlock [g_, s_, v_, twoN_,prec J : =Module [ 
{Omega, Alpha, M,p, den, r,num} , 

(M=2*Ceiling[twoN/2] ;p=PadeApproximant [Exp [z] ,{z,0,{M-l,M}}] ; 
den=Denominator [p] ;r=Roots [den==0 ,z] ; Alpha=Table [r [ [i , 2] ] ,{i , 1 ,M}] ; 
num=Numerator [p] ; hospital =num/D [den, z] ; 

Omega=SetPrecision [Table [hospital/. z->Alpha[[i]] ,{i,l,M,2}] ,prec+50] ; 
Alpha=SetPrecision [Table [Alpha[[i]] ,{i,l,M,2}] ,prec]) ; 

SetPrecision[-(2/v) Sum [Re [Omega [[i]] g/ . s->Alpha[ [i] ] /v] ,{i,l,M/2}] ,30] 

] 

In the above algorithm, g = <?(s), s = s, v = v, twoN=27V in Eq. (|22p . and prec = precision of calculation. Typical 
values are twoN = 8 and prec = 70. The algorithm, which is quite fast, returns the numerical value of G(v). It utilizes 
that fact that the partial fraction approximation made for e z is equal to a Pade approximant whose numerator is 
order 2N — 1 and whose denominator is order 2N. 

The algorithm first insures that M=twoN is even. It then constructs p, the Pade approximant whose numerator 
is a polynomial of order M — 1 and denominator a polynomial of order M. It finds r, the complex roots of the 
denominator, which are cti, the poles of Eq. ([22]) . Using L'Hospital's rule, it finds the residue u>i corresponding to the 
pole on. At this point, all of the mathematics is symbolic. It next finds every other pair of (a,i,uji) to the desired 
numerical accuracy; they come consecutively, i.e., c\\ = a.2, 0J\ — (I>2, a 3 — 6l±, u>3 = 0)4, etc. Finally, it takes the 
necessary sums, again to the desired numerical accuracy, but only over half of the interval i = 1, 3, . . . , N , by taking 
only the real part and multiplying by 2. 

If g(s), the input to the algorithm, is an analytic relation and v is a pure number (from the point of view of 
Mathematical 31/10 is a pure number, but 3.1 is not), then, for sufficiently high values of prec, you can achieve 
arbitrarily high accuracy. If we define the accuracy as 1 — G(v) n umericai/G(u)truc, numerical tests on many different 
functions shows that it goes to for large 2N as an inverse power law in 2N. On the other hand, if g(s) is obtained 
numerically, cither from having to find the Laplace transform /(s) and/or h(s) from numerical integration techniques, 
the overall accuracy of inversion is limited by the need to only use relatively small values of 2N — in the neighborhood 
of 6 — 12, limiting the overall accuracy to be in the neighborhood of 10 -5 , which fortunately is ample for most 
numerical work. Typically, numerical integration routines arc not accurate to better than ~ 10~ 5 , and one can not 
use ui's — which alternate in sign — that are larger than ~ 10 14 — 10 16 , which occur for relatively small values of 2N. Of 
course, this is not a limitation if g(s) is able to be expressed in closed form. 



APPENDIX B: MATHEMATICA NUMERICAL LAPLACE TRANSFORM ALGORITHM 



We note from Eq. (jTJJ) that the function that we must invert is 

g(s) = f(s)/h(s). (Bl) 

Although we can obtain an analytic solution for h(s) — true for LO as well as NLO — we must find f(s) by numerical 
means. From Eq. (fT2")> . we see that the Laplace transform of F(v) is given by the integral 

/•OO 

f{s) = / F(v)e- sv dv. (B2) 
Jo 

Because of the 00 at the upper limit of the integral, for numerical work it is useful to make the transformation 
v = — lnu, yielding 

h(s)= [ H(-]nu)u s - 1 du. (B3) 



The upper limit of the integral is now finite, leading to the stable Mathematica algorithm: 

NLaplaceTransf orm [H_, vv_, ss_] := NIntegrate [u~ (ss - 1)*H /.vv -> -Log[u], {u, 0, 1}, 
AccuracyGoal -> 15, MaxRecursion -> 15] 
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APPENDIX C: GLOBAL PARAMETERIZATION OF F^ p (x,Q 2 ) USING ZEUS STRUCTURE FUNCTION 

DATA 

Berger, Block and Tan [l5[ have parameterized the proton structure function F2 P (x, Q 2 ) as 



F™(x,Q 2 ) = (l-x) 



F P 



1 — Xp 



A{Q 2 )\n 



xp 1 — x 



x 1 — Xp 



B(Q 2 )\n 2 



xp 1 — x 

X 1 — Xp 



(CI) 



Here xp = 0.09 specifies the location in x of an approximate fixed point observed in the data where curves for different 
Q 2 cross. At that point, dFj p (x P , Q 2 )/d\nQ 2 « for all Q 2 ; F P = F] p (x P ,Q 2 ) = 0.41 is the common value of F 2 7P . 
The Q 2 dependence of F^ p (x,Q 2 ) is given in those fits by 



A{Q 2 ) = a + a l \nQ 2 ±a 2 ln 2 Q 2 , 
B(Q 2 ) = 6 + 6i lnQ 2 ±6 2 ln 2 Q 2 . 



(C2) 



The fitted quantities and their errors are shown in Table |TJ 



TABLE I: Results of a 6-parameter fit to ZEUS F p (x, Q 2 ) structure function data [T(| using the x and Q 2 behaviors of Eq. (|C1|) 
and Eq. (|C2|I . with Q 2 in GeV 2 . The renormalized Xmin P er degree of freedom, taking into account the effects of the Ax? max = 6 
cut [ljj, is given in the row labeled TZ x Xmin/d.f. The errors in the fitted parameters are multiplied by the appropriate r X 2 [13] ■ 



Parameters Values 



a 




-5.381 x 10" 2 ±2.17 x 10 


-3 


a± 




2.034 x 10~ 2 ± 1.19 x 10" 


3 


a 2 




4.999 x 10" 4 ± 2.23 x 10" 


4 


bo 




9.955 x 10 -3 ± 3.09 x 10" 


4 


bi 




3.810 x 10~ 3 ± 1.73 x 10" 


4 


62 




9.923 x 10~~ 4 ±2.85 x 10" 


5 


Xmin 




165.99 




K x 


Xmin 


184.2 




d.f. 




169 




TZ x 


Xmin/d.f 


1.09 



The fit to the data on i^ p (x, Q 2 ) was restricted to the region x < xp. In the absence of a global fit to the data in 
this region, they simply extended their parameterization, piecewise, to the large-x region, using the form 



F^(x,Q 2 ) = F P (f-J [T^f;) 



1 



where the piecewise extension on the r.h.s. of Eq. (|C3[) is obviously continuous with the l.h.s. at x = xp, independently 
of fi(Q 2 ). The exponent fi{Q 2 ) is determined by requiring that the first derivatives with respect to x of the function in 
Eqs. (|C1|) and the function on the r.h.s. of (|C3[) also match at x — xp. These results give the required parameterization 



of F2 P (x, Q 2 ) over all x, required in Eq. (JH to evaluate T(x, Q 2 ), and eventually, T(v, Q 2 ) of Eq. (JT)), needed to evaluate 
the exact solution for G(v, Q 2 ) in Eq. (JTHJ) , as well as calculating the Laplace transform /(s, Q 2 ) used in the numerical 
solution of Eq. (fi"4")) . which was calculated using the numerical algorithm of Appendix IBl 
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